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Temperature fluctuations in the normal direction of planar crystals such as graphene are quite 
violent and may be expected to influence strongly their melting properties. In particular, they 
will modify the Lindemann melting criterium. We calculate this modification in a self-consistent 
Born approximation. The result is applied to graphene and its wrapped version represented by 
single- walled carbon nanotubes (SWNTs). It is found that the out-of-plane fluctuations dominate 
over the in-plane fluctuations. This makes strong restrictions to possible Lindemann parameters. 
Astonishing we find that these large out-of-plane fluctuations have only a small influence upon the 
melting temperature. 

PACS numbers: 61.46.Fg, 64.70.D-, 68.60.Dv, 87.16.D- 
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The production of macroscopic two-dimensional (2D) 
graphene sheets by mechanical cleavaging [l[ has demon- 
strated that free-standing or suspended 2D crystals can 
exist despite large positional fluctuations in two dimen- 
sions. Since then, a variety of other free-standing 2D 
crystallites have been prepared 0. Wrapped versions 
of the 2D free-standing graphene had been found much 
earlier in 1991 [|[. Recent observations Q have confir- 
mend the theoretical expectation that freely suspended 
graphene sheets are strongly undulated and behave more 
like solid membranes than two dimensional (2D) crystals 
0. The undulations are a consequence of the thermal 
fluctuations of the membrane. In this paper, we calculate 
these fluctuations quantitatively and discuss their impli- 
cations upon the melting properties such as Lindemann 
parameter and melting temperature. Q The results will 
be compared with corresponding 2D crystals. 

The easiest way to estimate the melting temperature 
of a three dimensional (3D) crystal is based on the Linde- 
mann criterium [||. According to it, a 3D crystal starts 
to melt when the square root of the the rmal ex pectation 
value of lattice site elongations Oi = v(uJ)t exceeds a 
certain fraction of the nearest-neighbor lattice distance a, 
usually around 0.1 — 0.15 Q- Above the melting temper- 
ature, the shear modulus of the lattice vanishes leading 
to a divergence in the displacement fluctuations typical 
for the liquid state. 

In two dimensions (2D) this criterium is no longer ap- 
plicable since the displacement fluctuations are always 
logarithmically divergent, reflecting the fact that after 
a long time, a 2D crystal migrates through the entire 
2D-space. There exists, however, a simple modification 
Q. Instead of (u?)x one may use the finite cumulants 
(^ij)T — ( r ij)T) where ry is the difference vector between 
the atoms associated with the nearest-neighbor lattice 
sites i and j. This leads to a modified Lindemann num- 
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Here A/i denotes the set of all nearest-neighbor lattice 
pairs and |Ai | their number. For the Lennard-Jones and 
Wigner lattices, Bedanov et al. § found by computer 
simulations values of £ s { ~ 0.15 — 0.2. We have de- 
rived the same values analytically for a triangular gen- 
eralization B of a square lattice defect melting model 
0. 

At this point it is useful to realize that a migration 
problem and an associated divergence of Oi exists also 
in three dimensions if the system is finite, i.e., for 3D 
clusters and polymers. There one defines a modified Lin- 
demann number 
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Here A2 is the set of all lattice site pairs whose num- 
ber is N(N — l)/2 where N is the number of atoms in 
the lattice. The number L\ was introduced by Kaelberer 
and Etters [lj|, the number £| by Berry et al. [Till . 
For small clusters, C\ and L\ have similar values [12| of 
around L\ 2 w 0.03 — 0.05 at the melting point. Above 
the melting point all modified Lindemann numbers in- 
crease considerably, but do not go to infinity (this being 
in contrast to <Zj). 

The main difference between C S ^ D and £f comes from 
the last term in the square root in (JTJ and ([2|). Whereas 
{?ij)T is the temperature average of the difference vector 
of sites i and j, i.e. = (xij , ytj , Zij) , the expectation 
value (rij)r is the average value of the bonding length of 



sites i and j, i.e. 
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one would expect that the 3D version of (JT|), jC^ 
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and L\ ([2]) could be equally useful in determining the 
melting point. This is indeed the case for 3D crystals. 

In this paper we shall consider all three Lindcmann 
numbers C^ D , C\, and L\ as candidates for a melting 
criterion for solid membranes such as graphene lattices or 
SWNTs. It will turn out that for these L S { 3D is unsuit- 
able for calculating the melting temperature. The reason 
lies in the large out-of-plane fluctuations of the mem- 
brane vary little when crossing the melting point. These 
fluctuations cancel each other in |2J since (4') T ^ and 

ry « (4 +Vii) 1/2 + (1/2)4/(4 +Vii) 1/2 but not in 
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Freely suspended graphene sheets are always undu- 
lated and behave like a solid membrane [J ■ Nelson et al. 
@j EH have shown that in-plane fluctuations tend to sta- 
bilize a solid membrane such that a flat phase can exist in 
spite of its large 2D fluctuations. The melting tempera- 
ture of (5,5) SWNTs was determined by Zhang et al. [3 
within numerical simulation to be around T m ~ 5000 K, 



in agreement with experimental determinations 15j . The 
value of the Lindcmann number £§ was around L\ « 0.03 
at the onset of melting defined by the abrupt increase of 
£3. However, when using the region of strong increase of 
the internal energy they obtain a range £fj ~ 0-03 ~ 0.07 
from the onset of melting to its completion. 

The shapes of SWNTs near the melting tempera- 
ture are in general strongly deformed from a pure tube 
form. This leads to the conclusion that the 2D nearest- 
neighbor Lindemann number C S { 2D (fT|) is not a useful 
quantity for a melting criterion. One rather should use 
the Lindemann-like numbers @ or the 3D form of ([T]) 
which both respect the 3D rotational symmetry of the 
system. In the following, we shall first calculate L\ which 
for small clusters and small supcrcells agrees in molecular 
dynamic simulations with L\ . We shall restrict ourselves 
to the (5,5) SWNT so that we can compare our theoret- 
ical results with existing simulation data. We shall find 
that despite the large vertical fluctuations of the mem- 
branes, the Lindemann number ([2]) depends mostly on 
the in-plane fluctuations and provides us with a valuable 
melting criterium. This is not the case for the 3D- version 
of the Lindcmann number C\ having its reason in the fact 
that the out-of plane fluctuations are even larger than 
the in-plane fluctuations for SWNTs and graphene. Sur- 
prisingly, the melting temperature of SWNTs is modified 
only little by these large out-of-plane fluctuations at high 
temperatures. 
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and Uj are the lattice displacements in sy-plane, while 
/ is the out-of-plane displacement. The constant /i is 
the shear modulus, and A is the Lame constant. The last 
term in (U) with the constant kq accounts for the bending 
stiffness of the membrane. The quantity (Xij is an exter- 
nal stress source which will help us to calculate fl}, ([2]) 
from derivatives of the partition function with respect to 
<7jj . The line element on the membrane for small distor- 
tions is given by dl' 2 — dl 2 + 2uijdxidxj , where dl 
is the length of the undistorted planar surface, which we 
identify with the equilibrium lengths hj between sites i 
and j. Thus we calculate C\ 2 by first inserting the line 
element dl' in ([5]), and afterwards expanding the resulting 
expressions for small displacements. Thus we obtain 
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Here e y ' are the unit vectors pointing from site i to j. 
In deriving (f5|) we used a Taylor expansion of the elonga- 
tion differences between two lattice sites (gradient expan- 
sion). This is justified for the small elongation differences 
of neighboring atoms occurring in C{ and £f up to the 
melting regime. In the atom pairs summed over in £2, 
the approximation is good only for small clusters or for 
small supcrcells in molecular dynamic simulations. For 
infinite solid membranes, this is no longer the case. The 
contributions of the far separate pairs (i, j) cause a dras- 
tic decrease of the Lindemann parameter L\ with the size 
of the system. In a 2D crystal, for example, the widely 
separated pairs contribute terms which grow logarithmi- 
cally with the separation: lim/^^ (4')t — ( r ij)T ~ lniy. 
This eliminates £2 f° r determining the melting point. 

In order to calculate ([5]) we first integrate out the xy- 
lattice displacement fields Ui in the partition function, 
leading to an effective Hamiltonian H — Hh + H aG with 
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d 2 x^fi [(h f + K) 2 - hi] + i/coCViVJ)} , (6) 
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II. MEMBRANES 



The elastic energy of a solid elastic membrane in the 
flat phase is given by [l6| 
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and the energy densities 
hf = \pL(VifVmf), 
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In (O we have used the abbreviation jj, = + v) = E/2 
where E is the Young modulus and v = A/(2/z + A) the 
Poisson ratio. The calculation of the energy densities 
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(J?]) is simplified by the fact that only the transverse part 

P /m V '/ V ™/ With P /m = (&m ~ <ll<lm/q 2 ) of the OUt-of- 

plane fluctuations is relevant after integrating out the 
in-plane fields Mi Q. The transverse projections lead 
to a useful restriction of the relevant phase space when 
calculating Feynman diagrams. 



A. Self-consistent Born approximation 

We now treat the Hamiltonian ((6]) within the self- 
consistent Born-approximation (SCBA) corresponding to 
the Hartree-Fock approximation for the cigenfunctions. 
Other approximations to the Hamiltonian ^ have been 
used [171 ] to calculate the universal roughening exponents 
of the membrane, for example in Ref. [l8| an extension of 
SCBA. 

Denoting the inverse Green function of the /- 
fluctuations by G _1 (k) = K k 4 + £(fc), we obtain from 
© within the SCBA 
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where we take into account only the Fock-part of the 
SCBA. It was shown in Ref. [HI that the Hartree-terms 
do not contribute for free boundary conditions of the xy- 
clongations of the membrane. To do the integral in (|S]) 
we use a circular Brillouin zone k < fcsz, with fcsz = 
871"/ -\/3a 2 for the triangular Bravais lattice of SWNTs and 
graphene. The integral ((8]) can be carried out exactly for 
small k fl9j to obtain the first two terms in the expansion 
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where Ct is a temperature-dependent constant which 
turns out to be 
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The symbol T denotes the dimensionless temperature 
T = flk B T / (n k B z) 2 ■ 

The second coefficient is determined as follows. We 
assume that the truncated small- fc expansion ([9]) can be 
used for all k in the Brillouin zone, implying that the 
inverse Green function has the form 



(11) 



with K r = Ko + K s- We shall see below that this as- 
sumption is indeed justified. At low temperature where 
Ct <C ftr&BZ we determine n r by inserting (JTTJ) into ([8]) 
and evaluating the integral for E(fc) at the momentum 
k = Ct I K- r ■ This momentum regime is most relevant in 
the integrals over G which we have to calculate in the 
following in order to determine the generalized Lindc- 
mann parameters. Moreover, we will show below that 
flSJ) is then justified in good approximation for momenta 



even in the whole Brioullin zone. At higher temperatures 
where Ct S> «V^bz, we determine n r by integrating ((8]) 
at k = fcBZ- In both temperature regimes we carry the 
momentum integrations up to k = k B z, and obtain in 
either case a quadratic equation for n r , solved by 
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for Ct <C n r k B z, 
for C T > K r k BZ . 
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Our approximations are justified in Fig. 1 showing in 
the main plot the quantity G~ 1 (fc) = K r k 4 + Crk 3 di- 
vided by the sum of Kok 4 and the numerically integrated 
right-hand side of the self-energy function ([5]) . The num- 
bers on the curves are the various dimensionless temper- 
atures T. Observing that the values of these curves are 
almost constant and equal to unity confirms that G(k) 
of Eq. (|11[) indeed fulfills almost exactly the SCBA equa- 
tion ||5J). The inset of Fig. 1 shows k t /ko as a function 
of the dimensionless temperature T calculated either by 
(fT2j) , corresponding in the figure to the (green) solid and 
dashed curves, or by the determination of n r / n$ by nu- 
merical integration of the right-hand side of ([8]) ((blue) 
dashed-dotted curve). The kink in this curve corresponds 
to parameter values where Ct / K r k B z = 1- Note that for 
graphene and SWNTs we have T m w 1.34 at the melting 
point T m w 5000 K. 

Next we calculate the expectation value (Vj/Vj/) 
where the average is taken with respect to the Gibbs 
measure of the Hamiltonian or ©, respectively. In 
SCBA, this leads to 
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Recalling ^ we observe that the strain in the xy-planc 



(10) is on the average equal to the negative of (fT3|): (V 
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ViUj)r = — (Vi/Vj/)r, implying that the self-induced 
stress due to thermal out-of-planc fluctuations vanishes. 



B. Lindemann numbers 

We are no prepared to calculate (£J) 2 of Eq. §5§ by dif- 
ferentiating the partition function of the elastic Hamil- 
tonian H c \ twice with respect to the stress source <Jij, 
and setting Uy = at the end. Going over to the 
effective Hamiltonian ([6]) we obtain two contributions 



(A) = oc; 



C zl where the first is the square of 



the Lindemann number ([TJ) for the 2D hexagonal solid 
given by 
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we obtain for (1 15[) with (fTB"]) by integration [T| 
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FIG. 1: (Color online) Upper figure shows G _1 (fc) » C T fc 3 + 
K r fe 4 divided by the sum of Kok 4 and the numerical integrated 
right-hand side of the self-energy function (J5J) where we used 
G~ 1 (k) = Crk 3 + K r k 4 . The numbers at the curves are the di- 
mensionless temperatures T. The inset shows the T-behavior 
of ft r /«o- The low-temperature dashed curve represents 
Eq. (|12p in the regime Ct <C K- r ksz, the high-temperature 
solid curve in the regime Ct 3> n, r ksz- The dashed-dotted 
curve pictures the ratio Kt/kq obtained from (JH} by numerical 
integration of the right-hand side with G~ 1 (k) — Crk 3 + K r k 4 . 



The average {• • -)id is calculated with respect to the 
Gibbs measure of the Hamiltonian ([3]) with / = corre- 
sponding to the 2D crystal. This contribution to {CI) 2 
comes from the derivatc of H aa in Eq. ([6]) with respect 
to <Tij. The second contribution C 2 zl has its origin in the 
derivates of the hh term in J6J with respect to <xy, and 
is found to be 
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where I(q) is the Fourier transform of the correla- 
tor (h f (x)h f (x')) T - (h f (x)) T {h f (x')) T . Within self- 
consistent Born approximation we have 
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The vertex correction factor A(k+q, k) is required within 
the SCBA by charge-current conservation. We first cal- 
culate (|16|) in the lowest approximation A w 1, to be jus- 
tified below. By using the analytic approximation (lll|) 
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^l + i>) 2 = (l + ^) 2 {l-2[*V(l + ^-^/(l + ^) 2 ]}. 
Our results depend strongly on the number of the 2D 
Young modulus E = 2jl. The literature gives a broad 
range of possible /i-values (see for example [2(| Hl[ and 
references therein) which makes the comparison of our re- 
sults with experiment non-straightforward. It was shown 
by Hsieh et al. [2fj| using a molecular dynamics simula- 
tion that the Young modulus of (5,5) SWNTs is soft- 
ened near the melting temperature to around 70% of 
the T = -value, in agreement with Dereli et al. (2l| 
where a simulation of the larger (10,10) SWNT was car- 
ried out. This value for the temperature reduction is 
in accordance with the temperature reduction of the 2D 
Young modulus in Wigner crystals [22j at melting deter- 
mined by computer simulation. This can be generalized 
by theoretical arguments to softening expressions of clas- 
tic moduli 0, El for 2D crystals in general. The T = 
-value for the (5,5) SWNT determined by Hsieh et al. 
[2fj ] is E ~ 660N/m and lies at the upper end of exist- 
ing Young moduli in the literature. On the other hand 
the simulation of Dereli et al. for the (10,10) SWNT 
results in a much lower value at room temperature of 
around E ps lAON/rn where it should not much differ to 
its T = -value [2(| lying at the lower end of existing 
Young moduli for SWNTs in the literature. One should 
compare this value with the value E ps 440 N/m found 
by Hsieh [2(| for the (17,0) tube taking into account that 
the Young modulus depends only on the diameter of the 
tube and not the helicity [2fj ] in first approximation. The 
origin of these discrepancies in the Young modulus values 
shown in the literature in general is not clear yet. 

To compare our analytic results with the simulation 
results of Zhang et al. [l4| we shall use in the follow- 
ing the T = -value E ps 350 N/m for the (5,5) SWNT 
which is in the immediate proximity of several simula- 
tions (see [23| and references therein) and experimental 
values [24j | . The associated softened 2D Young modu- 
lus for (5,5) SWNTs is thus E ps 245N/m, which will 
be used in the rest of the paper. The remaining param- 
eters are less significant for our results. We shall take 
v ps 0.14, k ps 6.24- 1 10- 18 Nm, fc 2 z ps 2.46 • 10 20 /m 2 
which are typical for SWNTs and graphene. Inserting 
these material parameters we obtain for the melting tem- 
perature [3 T m = 5000 K a value C\ ps 5.6 • 10" 4 . 
The contribution {HJ), on the other hand, adds to this 
{C S { 2D ) 2 ps 8.9 • 10~ 3 at T m ps 5000 K so that the modi- 
fied Lindcmann number C\ is '. 
in the the range of values C\ s 
numerical simulation [l4j|- 

Our calculation shows that the abrupt increase of C\ 
is a meaningful criterium for the determination of the 
melting point. At the melting point, the in-plane shear 
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modulus (i will drop to zero, where according to Eqs. (1141) 
and (fT7|) . Cf{ 2D goes to infinity. 

Next we calculate the 3D form of (J]) (Cl' 3D ) 2 = 
(C S { 2D ) 2 + C 2 zl + C 2 z2 + |(Vi/Vi/)T where the last term 
is due to nearest-neighbor out-of-planc fluctuations given 
by (|13p . The first three terms measure in-plane fluctua- 
tions where C 2 2 is given by the momentum integral 
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The functions Ifj(q) and lfjij(q) denote the Fourier 
transforms of the expectation values 

[^ l f{x)\/ J f{x)h f {x')) T - {^if(x)^ j f(x)) T {h f {x')) T ] 



1 



and 
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respectively. The contribution C 2 2 is calculated in the 
same way as C 2 zl [lj|, yielding 
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C T > n r k B z- 
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Here we have used the abbreviation (1 + £) = 1 + 
(1 + i/) [-0.31 + 0.25i>/(1 + i/)]. Using the material 
parameters of SWNTs and graphene given above, we 
find that the main contribution to Cl' 3D comes from 
the out-of-plane fluctuations and is given by C^' 30 w 
((V;/VJ) T /2) 1/2 w 0.22 at the melting point T m w 
5000 K. Thus, we find that the out-of-plane fluctuations 
((Vi/Vi/)T/2) 1 / 2 are even larger than the dominant 
contribution to the in-plane fluctuations C*' 2D . By com- 
paring the temperature dependence (fT3|) with (fT4|) we 
obtain that this is even the case for smaller tempera- 
tures. Furthermore, we realize that in contrast to the 
Lindcmann number C\ , the abrupt increase of the Lindc- 
mann number C S -^ 3D gives no good signal for the melting 
point of a solid membrane. The reason is that the vanish- 
ing elastic shear modulus \i at melting contributes in two 
ways to the dominant fluctuation term (Vi/Vi/)T/2 (fT3|) 
but neither of them changes this value much at melting. 
First, the out-of-plane fluctuations depend on \x via K r j ' Kq 
and remains finite for fi — > 0, and second they depend 
pick up logarithmic dependence on fi from Ct/ko^bz- 

Consider now the higher-order vertex corrections col- 
lected in the factor A(k + q, k) in Eq. (TTBT) . First we 
note that for A = 1 wc obtain 2p,I(q) /kgT < 3/8 in the 



dominant integration regime of ((15)) near q ps Ct/k for 
Ct ^ K r fcez and near q m /cbz for Ct ^> K r fcez- The 
factor 3/8 comes mainly from the reduction of the phase 
space integral by the projections P T in the polarisator. 
We expect that the nth order in A contributes roughly 
with a factor [2fiI(q)/kBT] n to I(q) in the dominant inte- 
gration regime of (|15p due to the additional phase space 
projection terms P T . We have checked this explicitly by 
taking into account first-order corrections in the vertex 
A in (1161) . A similar suppression of higher-order vertex 



correction contributions occurs in l\- (q) and 7?-^ (q) . 



C. Melting temperature 

Let us finally discuss the impact of the large out-of- 
plane fluctuations upon the melting temperature. In 
Ref. we have calculated the melting temperature of a 
2D triangular lattice approximately from the intersection 
of high- and low-temperature expansion of the free ener- 
gies associated with the Hamiltonian ([3]) with zero ver- 
tical displacements /(x). The transition was caused by 
integer- valued defect gauge fields accounting for the plas- 
tic deformations of the crystal in the xy-plane. These are 
coupled minimally to the xy-displacement fields Uj(x). 
In that theory, the melting temperature T m was found to 
obey the equation 
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where vp = y/3a 2 /2 denotes the 2D-volume (area) of the 
fundamental cell. In SWNTs and graphene, this result 
is modified by a factor S e~ 2W , where 5* is a structure 
factor and e~ 2W a Debye- Waller-like factor caused by 
the out-of-plane fluctuations /(x) ^ 0. The honeycomb 
lattice of SWNTs contains two atoms per triangular fun- 
damental cell leading to a structure factor S = 1/2. To 



estimate the size of 



-2W 



we observe that in the defect 



melting model, the defect gauge field appears at a similar 
place as the vertical distortion Vj/Vj//2 in the Hamil- 
tonian Q. Thus one can immediately write down the 
Hamiltonian Hh of Eq. (|6|) coming from defects. This 
leads to the low-temperature expansion of the partition 
function. In the high-temperature expansion, there ex- 
ist a dual stress representation of the partition function 
0, Q. In both low- and high-temperature representa- 
tions, the coupling terms between the defect fields or 
the stress fields to the out-of-plane fluctuations /(x) are 
smaller than the pure defect and stress term by approx- 
imately a factor (47T 2 fi) 2 C zl and {2ir) 2 C 2 z2l respectively. 
When neglecting these small coupling terms we find that 
the partition function receives a sizable correction fac- 
tor only in the the low-temperature approximation due 
to the Fock energy of the Hamiltonian H^. The Hartrec 
energy is missing as a consequence of the open-boundary 
conditions on the membrane 131. From these considcra- 
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Vf 1 
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Using the parameters above for (5,5) SWNTs we ob- 
tain W w 0.06 at T w 5000 K. The factor e~ 2W gives 
thus only a small correction to the melting tempera- 
ture determined by (|2H1) . The explicit evaluation of that 
equation yields a melting temperature T m w 8000 K 
(VF w 0.075), somewhat larger than the melting tem- 
perature T m f=a 5000 K of Zhang et al. [3| obtained by 
numerical simulation. 



III. CONCLUSION 



tions even at low temperatures. Thus they may be ex- 
pected to have dramatic consequences for the Lindcmann 
numbers as well as the melting temperature of solid mem- 
branes in comparison to 2D crystals. Surprisingly, for 
the melting temperature this expectation was not con- 
firmed. The fluctuations was discussed by evaluating the 
3-D-vcrsion £ 1 ' 3 ' D of the Lindemann number JT]), orig- 
inally introduced to estimate the melting temperature 
of 2D solids, and the Lindcmann number CI defined in 
Eq. (|2|), originally introduced in cluster physics. We ob- 
served that a Lindcmann critcrium based on C\ is more 
reliable than that based on the former. The associated 
Lindemann number is dominated by in-plane fluctua- 
tions, in contrast to the former which is dominated by 
the large out-of-plane fluctuations. By calculating, in 
addition, the melting temperature from a simple defect 
model of melting for single- walled carbon nanotubes and 
graphenc (j20|) we observed in contrast to the expectation, 
that the melting temperature depends only very little on 
the large out-of-plane fluctuations. 



In this paper we have calculated the fluctuations of 
solid membranes like graphene and single-walled car- 
bon nanotubes with the help of the self-consistent Born- 
approximation. Our results show that the out-of plane 
fluctuations are much larger than the in-plane fluctua- 
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